In [ ]:

    

from astropy.io import ascii
import astropy.units as u
import astropy.coordinates as coord
import numpy as np
from numpy.polynomial.polynomial import polyval
import matplotlib.pyplot as pl
pl.style.use('apw-notebook')
%matplotlib inline
import h5py
from scipy.ndimage import gaussian_filter
from scipy.stats import scoreatpercentile


    

In [ ]:

    

ps1_filename = "/Users/adrian/projects/globber/data/ngc5897/PS1_stars_pv3_dered_sm.npy"
iso_filename = "/Users/adrian/projects/globber/data/ngc5897/dartmouth_iso_ps1.dat"


    

In [ ]:

    

cluster_c = coord.SkyCoord(ra=229.352*u.degree,
                           dec=-21.01*u.degree)
DM = 15.55

color = ('g', 'i')
mag = 'i'


    

In [ ]:

    

# TODO: swap this out for just using the XCov file instead
ps1 = np.load(ps1_filename)
ps1_c = coord.SkyCoord(ra=ps1['ra']*u.degree, dec=ps1['dec']*u.degree)

mask = (ps1['iErr'] < 0.1) # & (ps1_c.separation(cluster_c) > 1.*u.arcmin)
ps1 = ps1[mask]
ps1_c = ps1_c[mask]


    

In [ ]:

    

# read dartmoth isochrone
iso = ascii.read(iso_filename, header_start=8)


    

In [ ]:

    

idx = (ps1_c.separation(cluster_c) < 6*u.arcmin)


    

In [ ]:

    

pl.figure(figsize=(6,6))
pl.plot(ps1['ra'][~idx], ps1['dec'][~idx], ls='none', marker='.')
pl.plot(ps1['ra'][idx], ps1['dec'][idx], ls='none', marker='.', color='g')
pl.xlim(cluster_c.ra.degree+0.5, cluster_c.ra.degree-0.5)
pl.ylim(cluster_c.dec.degree-0.5, cluster_c.dec.degree+0.5)


    

In [ ]:

    

x0 = ps1['dered_{}'.format(color[0])]-ps1['dered_{}'.format(color[1])]
m0 = ps1['dered_{}'.format(mag)] 

fig,axes = pl.subplots(1,3,figsize=(10,6),sharex=True,sharey=True)

axes[0].plot(x0[~idx], m0[~idx],
             ls='none', marker=',', alpha=0.04)

axes[1].plot(x0[idx], m0[idx],
             ls='none', marker=',', alpha=1.)
axes[1].plot(iso['{}P1'.format(color[0])]-iso['{}P1'.format(color[1])], iso['{}P1'.format(mag)]+DM,
             ls='-', marker=None, alpha=0.5, lw=3)

axes[2].plot(x0[~idx], m0[~idx],
             ls='none', marker=',', alpha=0.04)
axes[2].plot(x0[idx], m0[idx],
             ls='none', marker=',', alpha=1.)

axes[0].set_xlim(-0.75,1.25)
axes[0].set_ylim(22, 13)


    

In [ ]:

    

ps1['ra'].min(), ps1['ra'].max()


    

In [ ]:

    

fig,axes = pl.subplots(1,3,figsize=(10,6),sharex=True,sharey=True)

_ix1 = ps1['ra'] < 221
print(_ix1.sum())
_ix2 = np.random.permutation(_ix1.sum())[:10000]
axes[0].plot(x0[_ix1][_ix2], i0[_ix1][_ix2],
             ls='none', marker='.', alpha=0.25)

_ix1 = ps1['ra'] > 239
print(_ix1.sum())
_ix2 = np.random.permutation(_ix1.sum())[:10000]
axes[1].plot(x0[_ix1][_ix2], i0[_ix1][_ix2],
             ls='none', marker='.', alpha=0.25)

# ----------------------------------------------------

_ix1 = ps1['ra'] < 221
_ix2 = np.random.permutation(_ix1.sum())[:5000]
axes[2].plot(x0[_ix1][_ix2], i0[_ix1][_ix2], color='k',
             ls='none', marker='.', alpha=0.25)

_ix1 = ps1['ra'] > 239
_ix2 = np.random.permutation(_ix1.sum())[:5000]
axes[2].plot(x0[_ix1][_ix2], i0[_ix1][_ix2], color='k',
             ls='none', marker='.', alpha=0.25)


axes[0].set_xlim(-0.75,1.25)
axes[0].set_ylim(22, 13)


    

In [ ]: